home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Digital / DFT.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  15.8 KB  |  542 lines
  1. (*  :Title:    Discrete Fourier Transforms  *)
  2.  
  3. (*  :Authors:    Wally McClure, Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    This file allows the computation of the forward and
  7.         inverse discrete Fourier transform (DFT).
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Digital`DFT`  *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*
  33.     :History:   start            July 30, 1989
  34.         redirection        October 2, 1989
  35.         made into package    April 18, 1990
  36.         made into release    May 25th, 1990
  37.  *)
  38.  
  39. (*  :Keywords:    discrete Fourier transform  *)
  40.  
  41. (*
  42.     :Source:    'Discrete-Time Signal Processing', Alan V. Oppenheim
  43.              and Ronald W. Schafer
  44.  *)
  45.  
  46. (*  :Warning:    *)
  47.  
  48. (*  :Mathematica Version:  1.2 or 2.0  *)
  49.  
  50. (*  
  51.     :Limitation: The DFT is calculated by first attempting to calculate the 
  52.              DTFT of a function using the MyDTFT rule base.  At the end 
  53.          of the MyDTFT rule base, a call to the z-transform is made.  
  54.          If  the  function  does  not have a transform specified in 
  55.          either domain, then a nasty expression is returned.    
  56.  *)
  57.  
  58. (*
  59.     :Discussion:  The DFT is based on the DTFT.  A signal is made finite
  60.           in extent by multiplying it by ( Step[n] - Step[n - N] ),
  61.           where N is the number of points at which to evaluate the
  62.           DFT.
  63.  *)
  64.  
  65. (*
  66.     :Functions:    DFTransform
  67.         InvDFTransform
  68.  *)
  69.  
  70.  
  71.  
  72. (*  B E G I N     P A C K A G E  *)
  73.  
  74. BeginPackage[ "SignalProcessing`Digital`DFT`",
  75.           "SignalProcessing`Digital`DTFT`",
  76.           "SignalProcessing`Digital`DSupport`",
  77.           "SignalProcessing`Digital`InvZTransform`",
  78.           "SignalProcessing`Digital`ZTransform`",
  79.           "SignalProcessing`Digital`ZSupport`",
  80.           "SignalProcessing`Support`TransSupport`",
  81.           "SignalProcessing`Support`ROC`",
  82.           "SignalProcessing`Support`SigProc`",
  83.           "SignalProcessing`Support`SupCode`",
  84.           If [ TrueQ[$VersionNumber >= 2.0],
  85.            "Algebra`SymbolicSum`",
  86.            "Algebra`GosperSum`" ] ]
  87.  
  88.  
  89. If [ TrueQ[ $VersionNumber >= 2.0 ],
  90.      Off[ General::spell ];
  91.      Off[ General::spell1 ] ];
  92.  
  93.  
  94. (*  U S A G E     I N F O R M A T I O N  *)
  95.  
  96. DFTData::usage =
  97.     "Data-tag for the symbolic DFT."
  98.  
  99. DFTransform::usage =
  100.     "DFTransform[function, NumberOfPoints, TimeVariables, \
  101.     FourierVariables, options] returns the DFT of function, \
  102.     where function is a function of TimeVariables defined from \
  103.     0 to NumberOfPoints - 1. \
  104.     Note that only the first two arguments \
  105.     are required and that DFTransform calls DTFTransform."
  106.  
  107. Finish::usage =
  108.     "Finish is a data-tag for the DFT."
  109.  
  110. InvDFTransform::usage =
  111.     "InvDFTransform[function, NumberOfPoints, FourierVariables, \
  112.     TimeVariables, options] takes the inverse DFT of function, \
  113.     where function is a function of FourierVariables defined from \
  114.     0 to NumberOfPoints - 1. \
  115.     Note that only the first two arguments \
  116.     are required and that InvDFTransform calls InvDTFTransform."
  117.  
  118. KVariables::usage =
  119.     "KVariables is a data-tag for the DFT."
  120.  
  121. Start::usage =
  122.     "Start is a data-tag for the DFT."
  123.  
  124. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  125.  
  126.  
  127. Begin[ "`Private`" ]
  128.  
  129.  
  130. (*  displayFixUp  *)
  131. displayFixUp[ trans_, op_ ] :=
  132.     Append[ fixUp[op, trans],
  133.         TransformLookup -> Replace[TransformLookup, op] ]
  134.  
  135. (*  validvarQ  *)
  136. validvarQ[ x_Symbol ] := True
  137. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  138. validvarQ[ x_ ] := False
  139.  
  140.  
  141. (*  B E G I N     D F T  *)
  142.  
  143. (*  Evaluation of DFT operator  *)
  144. Unprotect[DFT]
  145. DFT/: TheFunction[ DFT[N_, n_, k_] [f_] ] := DFTransform[f, N, n, k]
  146. Protect[DFT]
  147.  
  148. (*  Options for DFTransform  *)
  149. DFTransform/: Options[ DFTransform ] := Options[ DTFTransform ]
  150.  
  151. (*  Extension of TheFunction to return the transform function.  *)
  152. DFTData/: TheFunction[ DFTData[ trans_, var_ ] ] := trans
  153.  
  154. (*  Magnitude/Phase plot of a DFT object *)
  155. mppoptions = {    Domain -> Discrete, DomainScale -> Linear,
  156.         MagRangeScale -> Linear, PhaseRangeScale -> Degree,
  157.         PlotRange -> All }
  158.  
  159. DFTData/: MagPhasePlot[ DFTData[trans_, KVariables[k_Symbol],
  160.                 Start[st_], Finish[end_] ] ] :=
  161.     MagPhasePlot[ trans, {k, 0, end - st}, mppoptions ]
  162.  
  163. DFTData/: MagPhasePlot[ DFTData[trans_, KVariables[{k1_Symbol, k2_Symbol}],
  164.                 Start[{s1_, s2_}], Finish[{e1_, e2_}] ] ] :=
  165.     MagPhasePlot[ trans, {k1, 0, e1 - s1}, {k2, 0, e2 - s2}, mppoptions ]
  166.  
  167. (*  This Code decides between a one-dimensional or multidimensional DFT. *)
  168. DFTransform[ f_ ] :=
  169.     Message[ Transform::novariables, "N (length)", GetVariables[f] ]
  170.  
  171. DFTransform[ f1_, N_ ] :=
  172.     DFTransform[ f1, N, DummyVariables[Length[N], Global`n] ]
  173.  
  174. DFTransform[ f1_, N_, n_ ] :=
  175.     DFTransform[ f1, N, n, DummyVariables[Length[n], Global`k] ]
  176.  
  177. DFTransform[ fun_, N_, n_, k_, options___ ]:=
  178.     Message[ Transform::badvar, "discrete-time", DFTransform, n ] /;
  179.     ! validvarQ[n]
  180.  
  181. DFTransform[ fun_, N_, n_, k_, options___ ]:=
  182.     Message[ Transform::badvar, "discrete-frequency", DFTransform, k ] /;
  183.     ! validvarQ[k]
  184.  
  185. DFTransform[ f1_, N_, n_, k_, options___ ] :=
  186.     Block [    {begin, end, newfun, notenough, op, trans, vars, w, wvars},
  187.         notenough = ( Length[k] < Length[n] );
  188.         If [ notenough,
  189.              Message[Transform::notenough, "k (DFT index)"] ];
  190.         vars = If [ notenough, DummyVariables[Length[n], Global`k], k ];
  191.         wvars = DummyVariables[Length[n], w];
  192.  
  193.         op = ToList[options] ~Join~ Options[DFTransform];
  194.         trans = If [ Length[n] > 1,  
  195.                  multiDDFT[f1, N, n, vars, wvars, op],
  196.                  MyDFT[f1, N, n, vars, wvars, op] ];
  197.  
  198.         begin = Start[ If [AtomQ[n], 0, Table[0, {Length[n]}]] ];
  199.         end = Finish[ N - 1 ]; 
  200.  
  201.         DFTData[FourierSimplify[trans], begin, end, KVariables[vars]] ]
  202.  
  203. (*  multiDDFT --  multidimensional DFT  *)
  204. multiDDFT[ fun_, Num_, varin_, varout_, wvars_, op_ ] :=
  205.     Block [    {F2 = fun, jj, length, num, var, w},
  206.         length = Length[varin];
  207.         For [ jj = 1, jj <= length, jj++,
  208.               var = varin[[jj]];
  209.               num = Num[[jj]];
  210.               w = wvars[[jj]];
  211.               F2 = MyDFT[F2, num, var, varout[[j]], w, op] ];
  212.         F2 ]
  213.  
  214. (*  MyDFT  *)
  215. MyDFT[ x_, N_, n_, k_, w_, op_ ] :=
  216.     DualDFT[ x, N, n, k, w, op, DFTransform ]
  217.  
  218. (*  E N D     D F T  *)
  219.  
  220.  
  221. (*  B E G I N     I N V E R S E     D F T  *)
  222.  
  223. (*  Messages  *)
  224. InvDFTransform::badlength = "Conflicting lengths in inverse DFT: `` != ``."
  225.  
  226. (*  Evaluation of inverse DFT operator  *)
  227. Unprotect[InvDFT]
  228. InvDFT/: TheFunction[ InvDFT[N_, k_, n_] [f_] ] := InvDFTransform[f, N, k, n]
  229. Protect[InvDFT]
  230.  
  231. (*  Options for InvDFTransform  *)
  232. InvDFTransform/: Options[ InvDFTransform ] :=
  233.     displayFixUp[ InvDFTransform,
  234.               { Definition -> False } ~Join~ Options[InvDTFTransform] ]
  235.  
  236. (*  Choose whether to do the one or multi dimensional inverse DFT`s.  *)
  237. InvDFTransform[ DFTData[x_, Start[st_], Finish[end_], KVariables[k_]] ] :=
  238.     InvDFTransform[ x, end - st + 1, k ]
  239.  
  240. InvDFTransform[ f_ ] :=
  241.     Message[ Transform::novariables, "N (length)", GetVariables[f] ]
  242.  
  243. InvDFTransform[ DFTData[x_, Start[st_], Finish[end_], KVariables[k_]], N_ ] :=
  244.     If [ SameQ[ end - st + 1, N ],
  245.          InvDFTransform[ x, N, k ],
  246.          Message[ InvDFTransform::badlength, N, end - st + 1 ] ]
  247.  
  248. InvDFTransform[ f_, N_ ] :=
  249.     InvDFTransform[ f, N, DummyVariables[Length[N], Global`k] ]
  250.  
  251. InvDFTransform[ f_, N_, k_ ] :=
  252.     InvDFTransform[ f, N, k, DummyVariables[Length[k], Global`n] ]
  253.  
  254. InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
  255.     Message[ Transform::badvar, "discrete-frequency", InvDFTransform, k ] /;
  256.     ! validvarQ[k]
  257.  
  258. InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
  259.     Message[ Transform::badvar, "discrete-time", InvDFTransform, n ] /;
  260.     ! validvarQ[n]
  261.  
  262. InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
  263.     Block [    {notenough, nvars, op, w, wvars}, 
  264.         notenough = ( Length[n] < Length[k] );
  265.         If [ notenough, Message[Transform::notenough, "k (frequency)"]];
  266.         nvars = If [ notenough,
  267.                  DummyVariables[Length[n], Global`n],
  268.                  n ];
  269.         op = ToList[options] ~Join~ Options[DFTransform];
  270.         wvars = DummyVariables[Length[n], w];
  271.  
  272.         If [ Length[k] > 1, 
  273.              multiDInvDFT[ fun, N, k, nvars, wvars, op ],
  274.              MyInvDFT[ fun, N, k, nvars, wvars, op ] ] ]
  275.  
  276. (*   Multidimensional inverse DFT  *)
  277. multiDInvDFT[ fun_, N_, k_, n_, w_, op_ ] :=
  278.     Block [    {F1 = fun, jj, length},
  279.  
  280.         length = Length[Num];
  281.         For [ jj = 1, jj <= length, jj++,
  282.               F1 = MyInvDFT[ F1, N[[jj]], k[[jj]],
  283.                      n[[jj]], w[[jj]], op ] ];
  284.  
  285.         F1 ]
  286.  
  287. (*  MyInvDFT  *)
  288. MyInvDFT[ x_, N_, k_, n_, w_, op_ ] :=
  289.     DualDFT[ x, N, k, n, w, op, InvDFTransform ]
  290.  
  291. (*  E N D     I N V E R S E     D F T  *)
  292.  
  293.  
  294. (*  B E G I N     D U A L     D F T  *)
  295.  
  296. DualDFT[ x_, len_, n_, k_, w_, op_, flag_ ] :=
  297.     Block [    {newdualrules, newexpr, newx = x, oldexpr = Null, ret, trace},
  298.  
  299.         (* generate the DFT rules *)
  300.         newdualrules = TransformFixUp[ n, k, op, dualDFT, False,
  301.                            DFTransform, Null, Null ];
  302.         If [ SameQ[flag, InvDFTransform],
  303.              newdualrules = newdualrules /. k -> -k ];
  304.         DualRules = Join[ newdualrules, OriginalDualRules ];
  305.  
  306.         (* determine the level of dialogue *)
  307.         trace = SameQ[ Replace[Dialogue, op], All ];
  308.  
  309.         (* rewrite trig functions in terms of exponentials *)
  310.         (* and put exponentials in a cannonical form       *)
  311.         newx = x /. TrigToExpRules;
  312.         newx = newx /. ( Exp[a_] :> Exp[ Factor[a] ] );
  313.         If [ trace && ! SameQ[x, newx],
  314.              Print[x];
  315.              Print["which is rewritten in a more convenient exponential form as"];
  316.              Print[newx] ];
  317.  
  318.         (* check for inverse DFT *)
  319.         newexpr = dualDFT[newx, n, k, w, len, flag, fixUp[op, flag]];
  320.         If [ SameQ[flag, InvDFTransform],
  321.              If [ trace,
  322.               Print[ "Using the Duality Theorem, the inverse DFT ",
  323.                    "of length ", len, " is the DFT of length ",
  324.                    len, " divided by ", len, " with the index ",
  325.                    "variable ", k, " negated in the result:" ] ];
  326.               newexpr = Rev[k][newexpr / len] ];
  327.  
  328.         (* take the 1-D dual DFT *)
  329.         While [ ! SameQ[newexpr, oldexpr],
  330.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  331.             oldexpr = newexpr;
  332.             newexpr = oldexpr /. ( dualDFT[a__] :>
  333.                     Replace[dualDFT[a], DualRules] ) ];
  334.  
  335.         (* fix up the result *)
  336.         newexpr = posttransform[oldexpr];
  337.         While [ ! SameQ[newexpr, oldexpr],
  338.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  339.             oldexpr = newexpr;
  340.             newexpr = posttransform[oldexpr] ];
  341.  
  342.         newexpr ]
  343.  
  344. (*  circularShiftValue  *)
  345. n0checkvalue = Null
  346.  
  347. n0checkQ[n0_] :=
  348.     Block [    {},
  349.         If [ ValueQ[n0checkvalue], n0checkvalue = n0 ];
  350.         SameQ[n0, n0checkvalue] ]
  351.  
  352. circularShiftQ[x_, n_, L_] :=
  353.     Block [    {demodx, modtoken},
  354.         Unset[n0checkvalue];
  355.         demodx = x /. ( Mod[n + n0_, L] :> modtoken /; n0checkQ[n0] );
  356.         FreeQ[demodx, n] ]
  357.  
  358. circularShiftValue[x_, n_, L_] := n0checkvalue /; circularShiftQ[x, n, L]
  359.  
  360. (*  fixUp  *)
  361. fixUp[op_, DFTransform] :=
  362.     { Definition -> Replace[Definition, op],
  363.       Dialogue -> Replace[Dialogue, op] }
  364.  
  365. fixUp[op_, InvDFTransform] :=
  366.     { Definition -> Replace[Definition, op],
  367.       Dialogue -> Replace[Dialogue, op],
  368.       Terms -> Replace[Terms, op] }
  369.  
  370. (*  posttransform  *)
  371. posttransform[x_] := x /. postrules
  372.  
  373. postrules = {
  374.  
  375.     Rev[v_][x_] :>
  376.         If [ Count[{x}, f_[p___][a___], Infinity] > 0,
  377.              Rev[v][posttransform[x]],          (* operators present *)
  378.              posttransform[x] /. v -> -v ] ,    (* no operators *)
  379.  
  380.     Rev[k_][ Rev[k_] [ x_ ] ] :> x ,
  381.  
  382.     CircularShift[n0_, L_, n_][a_. Impulse[n_ + b_.]] :>
  383.         a Impulse[n + b + n0] /;
  384.         FreeQ[{a, b}, n] ,
  385.  
  386.     dualDFT[ x_, n_, k_, w_, L_, DFTransform, op_ ] :>
  387.         Block [    {dtft},
  388.             dtft = DTFTransform[ x Pulse[L, n], n, w, op ];
  389.             RewriteSummations[ TheFunction[dtft], w, 0, 2 Pi ] /.
  390.                 w -> (2 Pi k / L) ],
  391.  
  392.     dualDFT[ x_, k_, n_, w_, L_, InvDFTransform, op_ ] :>
  393.         InvDTFTransform[x /. (k -> w L / (2 Pi)), w, n, op] ,
  394.  
  395.     Impulse[-v_Symbol] :> Impulse[v]
  396. }
  397.  
  398.  
  399. Format[ dualDFT[ x_, n_, k_, w_, L_, rest___ ] ] :=
  400.     SequenceForm[ ColumnForm[ { "DFT",
  401.                  StringJoin["   ", ToString[L], ", ", ToString[n]] } ],
  402.               { x } ]
  403.  
  404. Format[ w ] := "w"
  405.  
  406.  
  407. DualRules = { }
  408.  
  409. OriginalDualRules = {
  410.  
  411. (*  Constant value--- some call it noise  *)
  412. dualDFT[ a_, n_, k_, w_, L_, flag_, op_ ] :> a L Impulse[k] /; FreeQ[a, n] ,
  413.  
  414. (*  Impulse  *)
  415. dualDFT[ a_. Impulse[n_ + n0_.], n_, k_, w_, L_, flag_, op_ ] :>
  416.     a Exp[-2 I Pi n0 k / L] /;
  417.     FreeQ[{a,n0}, n] && Implies[NumberQ[n0], IntegerQ[n0]] ,
  418.  
  419. (*  Multiplication by a complex exponential  *)
  420. dualDFT[ Exp[Complex[0, k01_] Pi k02_. (a_. n_ + b_.)] x_., n_, k_,
  421.         w_, L_, flag_, op_ ] :>
  422.     Block [    {cshift, dft},
  423.         cshift = a k01 k02 L / 2;
  424.         If [ ! IntegerQ[cshift] && SameQ[Replace[Dialogue, op], All],
  425.              Print[ "assuming ", cshift, " is an integer" ] ];
  426.         dft = dualDFT[ x, n, k, w, L, flag, op ];
  427.         Exp[-2 I Pi (b/a) k / L] CircularShift[cshift, L, k][dft] ] /;
  428.     FreeQ[{a, b, k01, k02}, n] ,
  429.  
  430. (*  Circular shift of sequence  *)
  431. dualDFT[ CircularShift[n0_, L_, n_][x_], n_, k_, w_, L_, flag_, op_ ] :>
  432.     Exp[2 Pi I k n0 / L] dualDFT[ x, n, k, w, L, flag, op ] /;
  433.     FreeQ[{n0, L}, n] ,
  434.  
  435. dualDFT[ x_, n_, k_, w_, L_, flag_, op_ ] :>
  436.     Block [    {circShift},
  437.         circShift = circularShiftValue[x, n, L];
  438.         newx = x /. ( Mod[n + circShift, L] -> n );
  439.         Exp[2 Pi I k circShift / L] *
  440.           dualDFT[ newx, n, k, w, L, flag, op ] ] /;
  441.     circularShiftQ[x, n, L] ,
  442.  
  443. (*  Conjugate of sequence  *)
  444. dualDFT[ Conjugate[x_], n_, k_, w_, L_, flag_, op_ ] :>
  445.     Conjugate[
  446.       CircularShift[0, N, k][
  447.         Rev[k][dualDFT[x, n, k, w, L, flag, op]] ] ] ,
  448.  
  449. (*  Reverse of a sequence  *)
  450. dualDFT[ Rev[n_][x_], n_, k_, w_, L_, flag_, op_ ] :>
  451.     Rev[k][ dualDFT[x, n, k, w, L, flag, op] ] ,
  452.  
  453. (*  Homogeneity  *)
  454. dualDFT[ a_ x_, n_, k_, rest__ ] :>
  455.     a dualDFT[ x, n, k, rest ] /;
  456.     FreeQ[a, n] ,
  457.  
  458. (*  Additivity  *)
  459. dualDFT[ x_ + y_, n_, k_, rest__ ] :>
  460.     dualDFT[ x, n, k, rest ] + dualDFT[ y, n, k, rest ] ,
  461.  
  462. (*  Apply the Definition  *)
  463. dualDFT[ x_, n_, k_, w_, L_, flag_, op_ ] :>
  464. Block [ {newx, state, sumhead, trans},
  465.     newop = fixUp[ Prepend[op, Definition -> False] ]; (* disable option *)
  466.     sumhead = If [ TrueQ[$VersionNumber >= 2.0], SymbolicSum, GosperSum ];
  467.     newx = ToDiscrete [ TheFunction[x] ];
  468.     trans = sumhead[newx Exp[-2 Pi n k / L], {n, 0, L-1}];
  469.     If [ SameQ[Head[trans], sumhead],
  470.          dualDFT[x, n, k, w, L, flag, newop],
  471.          TransformDialogue[Definition, op, sumhead, x, trans] ] ] /;
  472.     Replace[Definition, op]
  473.  
  474. }
  475.  
  476. (*  E N D     D U A L     D F T  *)
  477.  
  478.  
  479.  
  480.  
  481.  
  482. (*  T R I G     R E W R I T E     R U L E S  *)
  483.  
  484. (* Trigonometric to Complex Exponential Simplification rules *)
  485.  
  486. TrigToExpRules = {
  487.     Sin[ a__ ]    :>     ( Exp[ I a ] - Exp[ - I a ] ) / ( 2 I ),
  488.     1/Sin[ a__ ]    :>     1/( ( Exp[ I a ] - Exp[ - I a ] ) / ( 2 I ) ) ,
  489.     Cos[ a__ ]    :>    ( Exp[ I a ] + Exp[ - I a ] ) / 2,
  490.     1/Cos[ a__ ]    :>    1/( ( Exp[ I a ] + Exp[ - I a ] ) / 2 ),
  491.     Tan[ a__ ]    :>    Sin[ a ] / Cos[ a ],
  492.     Cot[ a__ ]    :>    Cos[ a ] / Sin[ a ],
  493.     Sec[ a__ ]    :>    1 / Cos[ a ],
  494.     Csc[ a__ ]    :>    1 / Sin[ a ]
  495. }
  496.  
  497.  
  498. (* Exponential rewrite rules *)
  499.  
  500. ExpToTrigRules = {
  501.     (a_. Exp[ Complex[0, b_] w_. ] + a_. Exp[ Complex[0, c_] w_. ] + x_.) :>
  502.         2 a Cos[Abs[b] w] + x /;
  503.         ( b == -c ),
  504.  
  505.     (a_. Exp[ Complex[0, b_] w_. ] - a_. Exp[ Complex[0, c_] w_. ] + x_.) :>
  506.         2 a Sin[b w] + x /;
  507.         ( b == -c ),
  508.  
  509.     (d_ - d_. Exp[ Complex[0, b_] c_. ]) :>
  510.         -2 I d Exp[ I b c / 2 ] Sin[ b c / 2 ],
  511.  
  512.     (d_. Exp[ Complex[0, b_] c_. ] - d_) :>
  513.         2 I d Exp[ I b c / 2 ] Sin[ b c / 2 ],
  514.  
  515.     (d_ + d_. Exp[ Complex[0, b_] c_. ]) :>
  516.         2 d Exp[ I b c / 2 ] Cos[ b c / 2 ]
  517. }
  518.  
  519. (*  E N D     T R I G     R E W R I T E     R U L E S  *)
  520.  
  521.  
  522. (*  E N D     P A C K A G E  *)
  523.  
  524. End[]
  525. EndPackage[]
  526.  
  527. If [ TrueQ[ $VersionNumber >= 2.0 ],
  528.      On[ General::spell ];
  529.      On[ General::spell1 ] ];
  530.  
  531.  
  532. (*  H E L P     I N F O R M A T I O N  *)
  533.  
  534. Combine [ SPfunctions, { DFTransform, InvDFTransform } ]
  535. Protect [ DFTransform, InvDFTransform ]
  536.  
  537.  
  538. (*  E N D I N G     M E S S A G E  *)
  539.  
  540. Print[ "The discrete Fourier transform (DFT) rule bases are loaded." ]
  541. Null
  542.